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Abstract 

Capturing the dynamics of granular flows at intermediate length scales can often be difficult. We 
propose studying the dynamics of contact networks as a new tool to study fracture at intermediate 
scales. Using experimental 3D flow fields with particle scale resolution, we calculate the time 
evolving broken-links network and find that a giant component of this network is formed as shear 
is applied to this system. We implement a model of link breakages where the probability of a link 
breaking is proportional to the average rate of longitudinal strain (elongation) in the direction of 
the edge and find that the model demonstrates qualitative agreement with the data when studying 
the onset of the giant component. We note, however, that the broken-links network formed in 
the model is less clustered than our empirical observations, indicating that the model refiects less 
localized breakage events and does not fully capture the dynamics of the granular fiow. 
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I. INTRODUCTION 



Forces in dense granular materials are transmitted through an intricate network of particle 
contacts. In jammed granular matter, this contact network is able to support stresses and 
strains up to a yield stress. However, under large enough stress or strain, the material will 
yield and deform, as particle contacts fail. This failure of the jammed state gives way to 
flow. The reproducible nature of many flow fields of granular matter indicates that 

the particle contacts fail as a continuum. However, it is not possible to fully leave the 
particle scale: shear band widths can be as small as five particle diameters, so unlike in solid 
mechanics or fluid dynamics, the scale of gradients in velocity is not well separated from 
the particle scale. Therefore continuum equations, for example the velocity field, cannot be 
fully justified. 

The dynamics of dense granular materials are also not amenable to a purely single particle 
view. Indeed, the physics of force networks as well as dynamics near jamming highlight the 
importance of the intermediate scale. However, actually measuring this intermediate scale 
is challenging: under shear strain, force correlations can be found for particles more than 
10 particle diameters apart sl, and cooperative structures such as force loops 4| have been 
found. The force correlations in the direction of the principal stress axis approximately 
follow a power law. String-like dynamics in which strings of particles are observed to move 
cooperatively do not exhibit a preferred length scale but rather a wide, power-law like 



distribution of lengths jsj. Elucidating collective dynamics in flows is further complicated 
since string-like dynamics are defined relative to a random fluctuating background, not 
relative to an overall flow. 

This intermediate scale is however ideal for observation through the lens of network 
theory, and we propose to utilize such an approach to study the onset of fracture in 3D 
structure. The heavy-tailed characteristics of both forces and string-like dynamics suggests 
no characteristic length scale exists. Instead of attempting to discover a characteristic length 
scale, our analysis produces the characteristic deformation of the whole system over which 
the contact network is substantially altered. 

Surprisingly little is known about the arrangement of grain contacts in a 3D jammed state 
in experiments, since grains are not transparent to standard imaging methods. Numerous 
studies have characterized the force network at the boundaries of 3D packings, for example 



Mueth et. al. used carbon paper to measure forces on the boundary of a cylindrical cell filled 
with glass beads and subjected to uniaxial compression Techniques for obtaining the in- 
terior contact network include freezing the arrangement with glue and cutting slices through 
the material or layer-by-layer removal of the grains, but these are destructive methods and 
are not suitable for investigations of system dynamics. 

Very recently, complex network methods have been used to analyze the computationally 
derived force network present in a granular material j7-9|. The focus of these studies is on 
the characterization of the instantaneous contacts and forces in the system. In this paper 
we introduce the time-evolving of the broken-links network as a new way to characterize 
particle rearrangement events in our experimental 3D granular flow. We utilize the language 
of percolation on networks to describe the dynamics of the network 10|, llll- This approach 
provides interesting results on the intermediate length and time scales that are difficult to 
capture with continuum or particle based approaches. The application of network theoretic 
concepts to this broken link network provides new insights at this intermediate length (and 
time) scale. Additionally, we implement a model of fracture based on the average rate of 
longitudinal strain and compare its results to our experimental data. 



II. 3D SHEAR FLOWS: IMAGING AND PROCESSING 

We determine the three dimensional microstructure of a jammed granular material using 



a laser sheet scanning approach we have described previously 12 1 which was adapted to 
study granular shear flows. The split-cell geometry as described in 13|] consists of a (15cm)^ 
cell with a 9cm diameter disk embedded in the bottom. The cell contains 5mm diameter 
acrylic beads immersed in an index matched triton and laser dye (Nile Blue 690 Perchlorate) 
solution, and is filled to a height of 10 particle diameters. The experiment is schematized in 
Figure [1^. 

Figure [T|d displays an example of an illuminated cross section of the granular fiow. We 
shear the system at a constant rate Q = 1.05 x 10~^ rad/s, pausing in three degree increments 
to recover the proximity network as a function of shear (and therefore time). By scanning 
a laser, we are able to image slightly more than half the shear cell and reconstruct the 
individual particle positions in 3D as a function of time. We bin particle velocities in both r 
(the radial distance from the center of the disk in cylindrical coordinates) and z (the height 



above the disk) and compute the average angular velocity u, and its derivatives with respect 
to r and z shown in Figures [TJ:; and[T]l. These quantities are used to develop a simple model 
for link breakage events due to longitudinal strain. The precision of the particle positions 
in the shear cell is approximately 2% of the radial cutoff defined below. 

In order to recover a network linking proximate particles, we need to define which particles 
interact as neighbors. However, it is not possible to detect with certainty whether particle 
pairs are touching - even a nm scale gap between particles would suffice to prevent force 
transmission between particles. Therefore we investigate, as a dynamic, statistical metric 
of neighbor interactions, how particles move with respect to nearby neighbors. Recent 



work 



14Nl6| has demonstrated that under small forcing, particles in a jammed system 



predominantly roll or slide past neighbors. The distribution in the angle, a, between the 
relative displacement vector of a pair of particles and the vector connecting their particle 
centers is peaked for tangential displacements, indicative of sliding or rolling behavior. For 
an ensemble of pairs of particles that do not touch and move independently the distribution is 
uniform, independent of cos (a) (for dynamics in three dimensions, for 2D motion it would be 
independent of a). Thus, by varying the threshold of radial displacements between particles 
and observing the distribution of the a cosines, P[cos(a)], we can determine an upper bound 
for particle contact. Figure [2] plots the peak value of the distribution of cos(a) vs. radial 
displacements. We note that there is a steep decline in the peak value of the distribution 
(and hence an increasing flatness in the distribution) for increasing radial displacement. 
By selecting the displacement immediately preceding the steep decline, we conservatively 
select a maximum threshold for particle contact. All particle pairs closer than the specified 
threshold then make up a proximity network. 



III. NETWORK ANALYSIS 

In order to capture the dynamics of the granular flow, we study the time-evolving broken- 
links network. A link is considered 'broken' relative to some reference frame if the particles it 
linked in the reference time frame move away from each other some time later (i.e. no longer 
meet the criterion for assigning contacts). Thus, the time evolving broken-links network is 
the collection of all such links (and the nodes/particles they are attached to). A schematic 
of this operation is displayed in Figure [3^. To lower the occurrence of false rearrangement 
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events, in order for a link to break, it must be broken in the current frame as well as the 
subsequent frame. We do not allow links to reform once they have been broken. 

We choose to focus on the broken network because we are primarily interested in the 
propagation of fracture events. However, in other applications, it may be instructive to 
consider the complimentary network of persistent links, defined as the set particles that are 
linked in the reference frame and remain in contact for all subsequent frames. Figure [3b 
illustrates the complementarity between the broken network and the persistent network. 

In addition to considering the experimentally determined broken network, we implement 
a model of link breakages where the probability of a link breaking is proportional to the 
average rate of longitudinal strain (elongation) in the direction of the edge. Using the actual 
experimentally measured positions of the particles, we begin by choosing a reference time 
and consider the contact network present at that time. To consider the effect of increasing 
shear, we take the experimental particle positions after the shear has been applied and for 
each edge in the reference network, find r' the position of the midpoint between the particle 
centers. Thus, r' serves as an approximation for the location of the edge in the shear cell, 
even if the particles have moved away from each other. 

Since the dominant average velocity is in the angular direction and uniform as a function 
of angle, we can approximate v = {vr,vg,Vz) = (0, ru{r, 2;), 0), and calculate both the strain 
rate tensor D, and the rate of longitudinal strain. 



6 = 6-1) -e (1) 

(2) 

P 

where f is taken relative to the center of the disk and e is the unit vector from one particle 
center to the other. For each amount of shear considered, we define the probability of 
breaking a link in the contact network, p, to be 

p oc e(e) (3) 

where is the Heaviside step function. Thus, we only assign probabilities of breaking to 
edges that have positive longitudinal strain (particles moving apart). Linear interpolation is 
used to evaluate the velocity gradients. The probability of breaking is normalized such that 
the expected value of the fraction of edges broken, f^, after each amount of shear is applied 
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is equal to the experimentally observed value, as seen in Figure S^. We also implement the 
model described above with noise. For a given amount of shear, a number Mb of the reference 
edges must be broken in order to match the experimentally observed fraction broken, ff,. We 
introduce a noise parameter < t] < 1 such that M^ri, of these edges are broken according 
to Eqini and the remaining (1 — r])Mi, edges are broken randomly with uniform probability. 

A particularly useful metric to capture the growth of the broken-links network is the size 
of its giant component. The giant component is the largest collection of connected nodes 
and edges present in the broken network at any given time. In network theory, it is the 



order parameter that signifies a phase transition in the behavior of the system [10|, The 
size of the giant component, Sg, is the fraction of nodes from the reference graph that are 
in the giant component. Thus, as shear is applied, one would expect the giant component 
of the broken network to grow in time. 

The optimal value of r] is found by maximizing the agreement between the model and 
the experimental data. A subset of the reference networks is used to calculate the giant 
component size as a function of shear, and r] is tuned to minimize the mean squared difference 
between the calculated and experimental curves. 

We consider one set of constant shear with a total of 240 frames, corresponding to 717 
degrees of rotation of the bottom plate. The first 210 frames are each taken as a reference 
frame and define 210 proximity networks. We allow particle tracks to be considered valid if 
the largest number of consecutive untracked frames for that particle is less than or equal to 
one. This condition is necessary since the experimental data contain some gaps, i.e. particles 
are not detected in some frames. For each network, the growth of the broken-links network 
and the giant component is studied over the subsequent 87 degrees, and we only consider 
particles whose tracks are valid for those 87 degrees. 

We average over all 210 reference networks. Figure H] plots the result of this analysis, 
with the size (number of nodes) of the giant component (relative to the reference network) 
as a function of the fraction of edges broken, ft,. It is clear that the fraction of edges broken, 
fb, plays the role of a tuning parameter, in the same way that an occupation probability 
is used to tune the behavior of a percolating system. The data suggest a continuous phase 
transition in the size of the giant component. One can map the fraction of edges broken 
to strain, via Figure H^. Thus, one can map the value of fb at which the giant component 
forms to a characteristic strain scale corresponding to the onset of global fracture, Schar ~ 15 
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degrees. 

Continuing the analogy with a percolation phase transition, we can look for a measure 
associated with a correlation length. To this end, we plot the average non-giant cluster size 
as a function of ft, in Figure |5^. Note that mean cluster size peaks approximately in the 
region in which the transition appears to occur in Figure |Dd. We can further investigate 
this transition by looking at the cumulative distribution of cluster sizes at ft values near 
this peak. The cumulative size distribution P{s' > s) is the probability that a cluster has 
a size s' greater than or equal to s. Figure \5jp plots the cumulative size distribution of 
clusters for three different ranges in fb. a fb range well before the transition, a range close 
to the transition (corresponding to values near the peak of (s)), and a range well after the 
transition. We immediately note that the distribution of cluster sizes for a range near the 
transition is much heavier tailed than both distributions before or after the transition, as 
one would expect for a percolation like transition. 

Additionally, we can investigate the structure of the broken-links network as it grows by 
computing the network's clustering coefficient, which reflects the extent to which neighbors 
of a node are themselves neighbors. The clustering coefficient for a node j is defined as 

_ mj) 

where kj is the degree of node j, the number of neighbors, and the number of triangles T{j) 
is the number of pairs of neighbors of j that are also linked to each other in the broken 
link network. The clustering coefficient measures the fraction of possible triangles centered 
at node j that are observed. The network clustering coefficient is simply the average over 
all nodes of the nodes' clustering coefficients C = {cj) Thus, a larger value of C for 

the broken-links network means that particle contacts tend to break more between nodes 
whose contacts have broken previously. Figure plots the average clustering coefficient as 
a function of the fraction of edges broken averaged over all runs. We note that both the 
model and the model with noise are less clustered than the data. Thus, while it appears 
that the model and model with noise appear to qualitatively reproduce the transition of the 
giant component, they fail to fully capture the structure of the data. 
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IV. DISCUSSION 



Shear zones and reproducible flow fields are a key feature of granular flows. In order to 
capture dynamics between particle and bulk scales, we introduced a new, complex networks 
approach to studying the dynamics of granular flow: the time-evolving broken-links network. 
A transition reminiscent of a phase transition in percolation is observed: a giant component 
of this network develops as shear is applied to this system. The appearance of this giant 
component reflects the onset of global fracture in the sense that the breaks in the original 
proximity network start to become globally connected. Additionally, we note a peak in 
the mean component size in proximity to the observed transition, as well as a heavy tailed 
distribution in the component size. This heavy-tailed length scale distribution is reminiscent 
of the one observed in string-like dynamics |5|. The network approach used here is ideal 
for studying granular phenomenon that exhibit a diverging length scale of this type. Our 
work shows that even in the case of a diverging length scale, it is possible to recover a 
characteristic strain scale corresponding to the onset of global fracture. We hypothesize that 
this strain scale can be linked to physical features of the system such as the loss of reversibility 
(this is the focus of ongoing work). Comparison to a model of link breakages, in which 
breakage events depend on the average longitudinal strain, highlights the collective nature 
of dense granular shear flows: the experimental data's broken-links network exhibits a larger 
clustering coefficient, indicating a more localized fracture than the model. This suggests 
that in a granular shear flow, breakage of a link between particles increases the likelihood 
of breakage of other nearby links. This further implies that the rearrangement events are 
cooperative. Indeed there is evidence for collective dynamics in dense granular flows, for 



example collective roll-like dynamics seen in chute flows 



18| and string like rearrangements 



in excited dense granular materials 19|. The network analysis provides a characteristic 
strain scale at which the network topology changes specifically without the need to focus on 
a particular lengthscale for collective dynamics. In conclusion, we believe that a dynamic 
networks approach offers a novel method to quantify granular flows at intermediate scales. 
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FIG. 1. (Color Online) (a) Schematic of our experimental setup, a 3D granular imaging method 
using a laser sheet to illuminate cross sections of the index matched, laser dyed shear cell fluid, (b) 
An example of an illuminated cross section of the granular flow, (c) and (d) ^ the r and 

z components of the numeric angular velocity gradient. The region of large gradients corresponds 
to an area close to the disk edge. Distances are normalized by the particle contact cutoff A* = 24.5 
(see Figure [2]). 
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FIG. 2. (Color Online) The peak value of the a cosines, -P[cos(a)] as a function of radial dis- 
placement Ar- Smaller peak values indicate more uniform distributions. Note the sudden decrease 
to uniformity as the radial displacement increases. We take the value after the sudden decrease 
(A* = 24.5) as our conservative upper bound on particle contact distance. Inset: The distribution 
of P[cos(a)] for touching A^ < 24.5 [circles] and non-touching A,. > 24.5 [squares] particles. 
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FIG. 3. (Color Online) (a) A schematic of the definition and growth of the broken-links network, 
b) Example cross section of the broken [green/light gray] and persistent [black] networks, plotted 
in real space. 
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FIG. 4. (Color Online) (a) The fraction of edges broken as a function of applied shear for the 
experimental data, averaged over all reference networks, (b) The size of the giant component Sg 
as a function of the fraction of broken edges for the experimental data [circles], longitudinal strain 
model [squares], and the longitudinal strain model with noise [t] = 0.82) [triangles]. Errorbars 
represent the standard error and are smaller than the markers. 
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FIG. 5. (Color Online) (a) The average non-giant cluster size as a function of the fraction of 
broken edges fb for the experimental data, (b) The cumulative size distribution for the non-giant 
components taken at three different ranges of 0.05 < < 0.06 [lower triangles], 0.11 < fb < 0.12 
[upper circles], 0.24 < fb < 0.25 [middle squares]. 
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FIG. 6. (Color Online) The clustering coefficient C, as a function of the fraction of edges broken, 
averaged over all reference networks. Bins are selected such that each bin contains the same number 
of data points. Error bars represent the standard error and are smaller than the markers. The 
experimental data [upper circles] are more clustered than both the model [squares] and the model 
with noise [triangles]. 
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